Motivation

The flow around an infinite, circular cylinder is a well-studied, but yet not completely understood problem in CFD (computational fluid dynamics). A simple exploration of the steady and unsteady flow at low Reynolds numbers is undertaken here, to provide a light overview of the problem.

Implementation

We implement all simulation with OpenFoam, analysis with Paraview and Python3, and documentation code in R Xie, Dervieux, and Riederer (2020).

Mesh Assembly : Preliminaries

We assemble a set of preliminary meshes for simulations with Reynolds number 20 and 110, consisting of an inner cylindrical segment and outer rectangular segment as follows. The following mesh schematics list dimension parameters and block configuration, courtesy of Dr.Fabrizio Bisetti, UT Austin.

Mesh Schema

The blocks and block vertices are as follows. The blue vertices denote the other Z-plane.

The edges are defined by listing points \(P\) as follows.

A symmetry plane is defined in OpenFoam s.t. boundary effects do not significantly affect results. So the resulting simulation will look as below.

For all the meshes unless otherwise noted, the parameters are as follows, with units in meters for all dimensions on the diagram, and Re denoting the Reynolds number.

  • f : blockMeshFactor (decreases cell dimensions in each direction by this factor), set by default : int(max(10, Re/3))
  • R : cylinder radius, default: 1/2
  • R2 : ring block outer radius, default: 3/2
  • H : height, default: 4
  • F : forward distance, default: 4
  • W : wake (backward) distance, default: 4 + Re*(1/15)
  • K : +/- distance in Z-axis — mostly irrelevant for anything in this project, default: 4

The preliminary meshes run_20_1 and run_110_1 simply follow the default settings.

Mesh Files

BlockMeshDict files are available here: run_20_1 and run_110_1, blockMesh logs are available here: run_20_1 and run_110_1,
checkMesh logs are available here: run_20_1 and run_110_1.

Mesh Figures

Mesh images are shown here, for detail comparision. We will show the run_20_1 mesh first and then the run_110_1 mesh.

  • run_20_1:

  • run_110_1:

Preliminary Solutions

We assemble a set of preliminary solutions for Reynolds number 20 and 110 (the previous meshes). The Re=20 mesh will be shown first.
Three contour plots will be shown: \(\frac{u}{U}\), \(\frac{v}{U}\), \(\frac{p}{\rho U^2}\), along with a streamline plot.
For Re=110 (the unsteady case), a time history of \(\frac{u}{U}\), \(\frac{v}{U}\), \(\frac{p}{\rho U^2}\) at \((x,y)=(5.5,\pm0.5)\) versus normalized time \(\frac{t}{\frac{D}{U}}\).

Note data for a couple of points (about 0.1%) was imputed due to file errors.

  • run_20_1:

  • run_110_1:

Time History plots for run_110_1:

Mesh Improvement

We can see that the mesh extent is fairly adequate, except possibly the height (H) for Re=110. All boundary flow is not perturbed, with the exception of the aforementioned region. We can also see vortices for the Re=110 case, and the inner block does not cause any boundary issues, so we will make only the following minor changes to each mesh.

For the mesh run_20_1, we will increase the blockMeshFactor from 10 to 15 (so a block increase of ~2.25). This mesh will be called run_20_2.

For the mesh run_110_1, we will increase height (H) from 4 to 6. This mesh will be called run_110_2.

Mesh Files

Files for the new meshes are listed below.

BlockMeshDict files are available here: run_20_2 and run_110_2, blockMesh logs are available here: run_20_2 and run_110_2,
checkMesh logs are available here: run_20_2 and run_110_2.

Mesh Figures

Mesh images are shown here, for detail comparision. We will show the run_20_1 mesh first and then the run_110_1 mesh.

  • run_20_2:

  • run_110_2:

Solutions

The Re=20 mesh will be shown first.
The same three contour plots will be shown: \(\frac{u}{U}\), \(\frac{v}{U}\), \(\frac{p}{\rho U^2}\), along with a streamline plot.
For Re=110 (the unsteady case), a time history of \(\frac{u}{U}\), \(\frac{v}{U}\), \(\frac{p}{\rho U^2}\) at \((x,y)=(5.5,\pm0.5)\) versus normalized time \(\frac{t}{\frac{D}{U}}\).

  • run_20_2:

  • run_110_2:

Time History plots for run_110_2:

Unsteady Flow and Strouhal Number

We now use the simulations run_110_3, run_110_4, and run_110_2 to calculate an approximate Strouhal number: \(St=\frac{f}{U/D}={fcn}(Re)\).

The sole difference between these runs is the blockMeshFactor parameter, which increases block density and decreases timestep.

Probes were placed at \((x,y)=(5.5,\pm0.5)\), and a Fourier Transform of the nondimensionalized y-velocity was taken with respect to nondimensionalized time, s.t. the highest amplitude resultant frequency is the Strouhal number. Agreement between Strouhal numbers calculated for each probe was verified.

A table of results follows:

  • Mesh | St | timestep
  • 110_3 0.180995 0.000666667;
  • 110_2 0.182648 0.000555556;
  • 110_4 0.181818 0.000444444;

So clearly the approximations converge to a Strouhal number around 0.182.

Plots for a probe from each mesh follow, along with the Fourier Transform.

  • Mesh 110_3

  • Mesh 110_2

  • Mesh 110_4

Looking at separation on Mesh 110_2, we see that it occurs at about 117 degrees, as predicted by Park et. al. (Numerical Solutions of Flow Past a Circular Cylinder at Reynolds Numbers up to 160, 1203).

The referenced plot from Park et. al.:

Steady Flow:

Radial and Tangential Velocity at Cylinder Wall

We compute radial and tangential velocities along lines \(AB\) normal to the cylinder wall, for angles \(\Theta = \frac{\pi}{4}, \frac{\pi}{2}, \frac{3\pi}{4}\), as seen in the diagram below. This is done for three meshes with identical dimensions but an increasing number of cells {8500, 19050, 34000, 53000}.

So for \(\Theta = \frac{\pi}{4}\), with increasing numbers of cells:




For \(\Theta = \frac{\pi}{2}\), with increasing numbers of cells:




For \(\Theta = \frac{3\pi}{4}\), with increasing numbers of cells:




A table of values comparing the relative recirculation length (L/D) with the components of the rate of strain tensor \(e_{rr},e_{r\theta}\) is shown below.
We define those components as follows.

  • (L/D) | \(e_{rr}\) at \(\pi/4\) | \(e_{rr}\) at \(\pi/2\) | \(e_{rr}\) at \(3\pi/4\) | \(e_{r\theta}\) at \(\pi/4\) | \(e_{r\theta}\) at \(\pi/2\) | \(e_{r\theta}\) at \(3\pi/4\) | Meshes in order from run_20_1 to 4.
  • 0.79809 | 0.1837 | 2.5977 | -0.4306 | 0.0586 | 1.4517 | 0.0965
  • 0.86399 | 0.1234 | 2.0425 | -0.3213 | 0.0399 | 1.8885 | 0.0778
  • 0.83199 | 0.0924 | 1.6610 | -0.2552 | 0.0299 | 2.1238 | 0.0651
  • 0.81322 | 0.0737 | 1.3937 | -0.2104 | 0.0239 | 2.2721 | 0.0553

We see that as the mesh resolution becomes finer, the values for \(e_{rr}\) and \(e_r\theta\) do converge, but not nearly as quickly as the velocity field does (as seen in the above plots)!

Note the L/D recirculation length is expected from the results in Park et. al. As seen in the below image, Park et. al. also predict a recirculation length of about 0.8 diameters.

Drag Coefficient for Re=20:

Calculation of force on the cylinder can be done by implementing a simple control volume approach, by measuring the far-field pressure and velocity at the inlet and outlet as follows. Note all units are non-dimensionalized, and the force is per unit length (of the cylinder).

By momentum conservation, we have the following, where y is the vertical dimension perpendicular to the cylinder axis and the fluid’s general velocity. Note the inlet-outlet notation denotes a pair of integrals for the inlet and outlet, subtracted from each other.

\[F=\int_{inlet-outlet} p \,dy + \int_{inlet-outlet} \frac{\partial{\vec{p}}}{\partial{t}} \]

\[\implies F=\int_{inlet-outlet} p \,dy + \int_{inlet-outlet} \rho\|{\vec{u}}\|\vec{u}_x \partial{y}\]

Note this is trivial to calculate, since all the values are non-dimensionalized by OpenFoam.

The resulting force is 1.358 in units of \(\rho U^{2}L\), which becomes a Cd of 2.717.

Comparing this to results from Tritton (Physical Fluid Dynamics, 1977), and Heddleson et. al. (Summary of Drag Coefficients of Various Shaped Cylinders), we expect a Cd around 2.4 for a Reynolds number of 20, and our result becomes reasonable. Reference images from Tritton and Heddleson are attached below.

References

Allaire J, Xie Y, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2021). Rmarkdown: Dynamic Documents for r. Retrieved from https://CRAN.R-project.org/package=rmarkdown
Bengtsson H (2021). R.utils: Various Programming Utilities. Retrieved from https://github.com/HenrikBengtsson/R.utils
Urbanek S (2021). Jpeg: Read and Write JPEG Images. Retrieved from http://www.rforge.net/jpeg/
Xie Y (2013). animation: An R Package for Creating Animations and Demonstrating Statistical Methods.” Journal of Statistical Software, 53(1), 1–27. Retrieved from https://doi.org/10.18637/jss.v053.i01
Xie Y (2014). “Knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing reproducible computational research, eds. V Stodden, F Leisch, and RD Peng,. Chapman; Hall/CRC. Retrieved from http://www.crcpress.com/product/isbn/9781466561595
Xie Y (2015). Dynamic Documents with R and Knitr 2nd ed. Chapman; Hall/CRC, Boca Raton, Florida. Retrieved from https://yihui.org/knitr/
Xie Y (2021a). Animation: A Gallery of Animations in Statistics and Utilities to Create Animations. Retrieved from https://yihui.org/animation/
Xie Y (2021b). Knitr: A General-Purpose Package for Dynamic Report Generation in r. Retrieved from https://yihui.org/knitr/
Xie Y, Allaire JJ, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman; Hall/CRC, Boca Raton, Florida. Retrieved from https://bookdown.org/yihui/rmarkdown
Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook. Chapman; Hall/CRC, Boca Raton, Florida. Retrieved from https://bookdown.org/yihui/rmarkdown-cookbook

Appendix

Thank you so much for reading this work!